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1.  Introduction 


Consider  the  NxN  system  of  llnesr  equations 

(1-)^  M X - b, 

where  the  coefficient  matrix  H Is  large,  sparse,  symmetric,  and  positive 
definite.  Such  systems  arise  frequently  In  scientific  computation, 
e.g..  In  finite  difference  and  finite  element  approximations  to  elliptic 
boundary  value  problems.  ^ this  report^wm  present  a package  of 
efficient,  reliable,  well-documented,  and  portable  FORTRAN  subroutines 
for  solving  these  systems. 

Direct  methods  for  solving  (1)  are  generally  variations  of 
symmetric  Gaussian  elimination.  We  form  the  U^DU  decomposition  of  A, 
where  U Is  unit  upper  triangular  and  D positive  diagonal,  and  then 
successively  solve  the  triangular  systems 

(2)  U^y“b,  Dz-y,  Ux^z. 

When  M Is  large  (N  » 1),  (dense)  Gaussian  elimination  Is  prohibitively 

3 

expensive  In  terms  of  both  the  work  ('  1/3  N multiplies)  and  storage 
2 

(N  words)  required.  But,  since  M Is  sparse,  most  entries  of  M and  U 
are  zero  and  there  are  significant  advantages  to  factoring  M without 
storing  or  operating  on  the  zeroes  appearing  In  M and  U.  Recently,  a 
number  bf  implementations  of  sparse  Gaussian  elimination  have  appeared 
based  on  this  idea,  cf.  [2, 5,6, 7]. 

In  section  2,  we  describe  the  scheme  used  for  storing  sparse 
matrices,  while  In  Section  3 we  give  an  overview  of  the  package  from  the 
point  of  view  of  the  user;  for  further  details  of  the  algorithms 
employed,  see  [4,9].  In  section  4,  we  Illustrate  the  performance  of  the 
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since  the  coefficient  matrix  M and  the  upper  triangular  factor  U 
are  large  and  sparse.  It  Is  Inefficient  to  store  them  as  dense  matrices. 
Instead,  we  store  matrices  using  a row-by-row  storage  scheme  used  In 
previous  Implementations  of  sparse  symmetric  Gaussian  elimination, 
cf.  [1,41. 


This  scheme  requires  three  one-dimensional  arrays:  lA,  JA,  and  A. 
The  nonzero  entries  of  M are  stored  row-by-row  In  the  REAL  array  A.  To 
Identify  the  Individual  nonzero  entries  In  a row,  we  need  to  know  In 
which  column  each  entry  lies.  The  INTEGER  array  JA  contains  the  column 
Indices  which  correspond  to  the  nonzero  entries  of  M,  i.e..  If  A(K)  - 
M(I,J),  then  JA(K)  - J.  In  addition,  we  need  to  know  where  each  row 
starts  and  how  long  It  Is.  The  INTEGER  array  lA  contains  the  Indices  In 
JA  and  A where  each  row  of  M begins,  i.e..  If  M(I,J)  is  the  first 
(leftmost)  entry  of  the  1^^  row  and  A(K)  “ M(I,J),  then  IA(I)  “ K. 
Moreover,  IA(N+1 ) Is  defined  as  the  index  In  JA  and  A of  the  first 
location  following  the  last  element  In  the  last  row.  Thus,  the  number 
of  entries  in  the  I*^^  row  Is  given  by  IA(I+1)  - IA(I),  and  the 
nonzero  entries  of  the  I^^  row  are  stored  consecutively  in 
A(IA(I)),  A(IA(I)+1),  ...,  A(IA(I+1)-1) 


k 


and  the  corresponding  column  Indices  are  stored  consecutively  In 


JA{IA(I)),  JA(IA(I)+l),  ...,  JA(1A(I+1)-1). 


1 


□ □ 
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For  exMple,  the  5x5  matrix 


1 0 2 3 0 

0 4 0 0 0 

2 0 5 6 0 

3 0 6 7 8 

0 0 0 8 9 


M - 


Is  stored  as 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

lA 

1 

4 

5 

8 

12 

13 

JA 

1 

3 

4 

2 

1 

3 

4 

1 

3 

4 

5 

4 

5 

A 

1 

2 

3 

4 

2 

5 

6 

3 

6 

7 

8 

8 

9 

Since  the  matrix  M is  symmetric,  it  suffices  to  store  only  the 
nonzero  entries  in  the  diagonal  and  strict  upper  triangle  of  M.  The 
storage  scheme  used  is  the  same  as  before  (except  that  nonzero  entries 
in  the  strict  lower  triangle  of  M are  Ignored),  and  our  example  matrix 
would  be  stored  as 


123456789 
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to  the  number  of  nonzero  entries  In  (the  dlsgonal  siul  strict  upper 
triangle  of)  M. 
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3.  A Sparse  Syetrlc  Matrix  Packaae 

The  package  consists  of  three  drivers  and  five  subroutines  (see 
Figure  1).  The  ordering  driver  (subroutine  ODRV)  uy  be  used  to 
symmetrically  reorder  the  variables  and  equations  so  as  to  reduce  the 
total  work  (l.e.  the  number  of  multiplies)  and  storage  required.  The 
solution  driver  (subroutine  SDRV)  Is  used  to  solve  the  (permuted)  system 
of  linear  equations.  The  test  driver  (program  STST)  sets  up  a model 
sparse  symmetric  positive  definite  system  of  linear  equations,  calls 
ODRV  to  reorder  the  variables  and  equations,  and  calls  SDRV  to  solve  the 
linear  system.  In  the  remainder  of  this  section,  we  describe  each  of 
these  routines  In  somewhat  greater  detail.  The  codes  themselves  are 
extensively  documented;  for  further  details  about  the  algorithms 
employed,  see  [4,9]. 

Figure  1:  A schematic  overview  of  the  sparse  symmetric  matrix  package 
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A.  The  Ordering  Driver  (ODRV) 

The  work  and  storage  required  to  solve  a large  sparse  systea  of 
linear  equations  clearly  depend  upon  the  zero-nonzero  structure  of  the 
coefficient  matrix.  But  since  this  matrix  Is  symmetric  and  positive 
definite,  we  could  equally  well  solve  the  permuted  system 

(3)  Q M y - Q b.  Q X - y 

given  any  permutation  matrix  Q.  The  permuted  system  corresponds  to 
symmetrically  reordering  the  variables  and  equations  of  the  original 
system,  and  the  net  result  can  often  be  a significant  reduction  In  the 
work  and  storage  required  to  form  the  U^DU  decomposition  of  A [3]. 

The  ordering  driver  (subroutine  ODRV)  uses  the  Important  heuristic, 
the  minimum  degree  algorithm  (implemented  in  subroutine  ORDER),  to 
select  Q.  ORDER  does  a symbolic  elimination  on  the  nonzero  structure  of 
the  system.  At  each  step,  it  chooses  a pivot  element  from  among  those 
unellmlnated  diagonal  matrix  entries  which  require  the  fewest  arithmetic 
operations  to  eliminate  (ties  are  broken  arbitrarily).  This  has  the 
effect  of  locally  optimizing  the  elimination  process  with  respect  to  the 
number  of  arithmetic  operations  performed.  See  [8,9]  for  more  details. 

ORDER  returns  two  one-dimensional  INTEGER  arrays  of  length  N:  P 
contains  the  permutation  of  the  row  and  column  indices  of  M,  i.e.,  the 
sequence  of  pivots;  and  IP  contains  the  inverse  permutation,  i.e., 
IP(P(I))  “ I for  I ••  1,2,...,N.  If  only  the  upper  triangle  of  M is 
being  stored,  then  the  representation  of  M (i.e.,  the  arrays  lA,  JA,  and 
A)  must  be  rearranged  using  the  subroutine  SRO  (PATH-2  in  ODRV). 
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The  user  may  bypass  ODRV  entirely  by  setting  P(l)  - IP(I)  - 1 
for  I - l,2.3,...,N.  Alternately,  the  user  may  substitute  another 
ordering  subroutine  for  ORDER,  as  long  as  It  produces  the  two 
permutations  P and  IP.  But  again.  If  only  the  upper  triangle  of  H Is 
being  stored,  the  representation  of  M must  be  rearranged  using  SRO. 

B.  The  Solution  Driver  (SDRV) 

The  solution  driver  (subroutine  SDRV)  is  used  to  solve  the 
(permuted)  linear  system.  Follo«#lng  Chang  [1],  SDRV  breaks  the  solution 
process  into  three  steps:  symbolic  factorization  (subroutine  SSF) , 
numerical  factorization  (subroutine  SNF) , and  back-solution  (subroutine 
SNS).  First,  SSF  determines  the  nonzero  structure  of  the  rows  of  U from 
the  nonzero  structure  of  the  rows  of  M.  Second,  SNF  uses  the  structure 
information  generated  by  SSF  to  compute  the  U*’DU  factorization  of  M. 
Third,  SNS  computes  the  solution  x from  the  factorization  generated  by 
SNF  and  the  right-hand  side  b. 

By  splitting  up  the  computation,  we  have  gained  flexibility.  To 
solve  a single  system  of  equations.  It  suffices  to  use  SSF,  SNF,  and  SNS 
(PATH-1  in  SDRV).  To  solve  several  systems  In  which  the  coefficient 
matrices  have  the  same  nonzero  structure,  it  suffices  to  use  SSF  only 
once  and  then  to  use  SNF  and  SNS  for  each  system  (PATH-2).  To  solve 
several  systems  with  the  same  coefficient  matrix  but  different  right 
hand  sides,  it  suffices  to  use  SSF  and  SNF  only  once  and  then  use  SNS 
for  each  system  (PATH-3). 
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C.  The  Test  Driver  (STST) 

The  test  driver  (proKram  STST)  Is  used  to  verify  the  perfonsance  of 
the  package  on  a particular  computer  system,  and  may  be  used  as  a guide 
to  understanding  how  to  use  the  package.  It  generates  the  coefficient 
matrix  for  the  standard  five-point  finite  difference  approximation  on  a 
3x3  grid  to  the  Poisson  equation  and  chouses  the  right-hand  side  so  that 
the  solution  vector  x Is  ( 1 , 2. 3, 4, S, 6, 7, 8, 9) . Since  H is  symmetric,  we 
can  specify  either  the  entire  matrix  (CASE>1)  or  only  the  upper  triangle 
(CASE-2).  STST  then  calls  ODRV  to  reorder  the  variables  and  equations, 
and  SDRV  to  solve  the  linear  system.  At  each  stage,  the  values  of  all 
relevant  variables  are  printed  out.  A sample  output  appears  as  Appendix 
4. 


T 
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4.  Perforaance 

One  of  the  moat  Important  aspects  of  any  package  Is  Its  performance 
In  terms  of  both  the  time  and  storage  required  to  solve  a typical 
problem.  In  Tables  1 and  II,  we  present  the  time  and  storage  required 
to  solve  the  familiar  five-point  finite  difference  equations  on  an  n x n 
grid  for  several  values  of  n.  These  computations  were  performed  In 
single  precision  on  an  IBM  370/158  using  the  FORTRAN  X optimizing 
compiler . 

Table  I 
Time  Required 

Five-Point  Operator  on  an  n x n Grid 


Grid 

1 STST 

[ ORDER 

SRO 

SSF 

SNF 

sec/* 

SNS 

Total 

20 

0.083 

0.657 

0.043 

0.100 

0.407 

14.016 

0.087 

0.593 

30 

0. 190 

1.750 

0.100 

0.257 

1.430 

13.002 

0.230 

1.917 

40 

0.3A0 

j 

3.656 

0.177 

0.503 

3.700 

J 

12.449 

0.470 

4.673 
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Table  II 
Storage  Required 

Five-Point  Operator  on  an  n x n Grid 


Grid 

A.  JA 

U 

1 JU 

L_ 

Total 

Multa. 

20 

1,160 

3,368 

1,889 

7,259 

35,195 

30 

2,640 

9*456 

4,538 

18,496 

127,666 

1 

40 

4,720 

19,926  ! 

8,423 

36,351 

334,937 

1 
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Sparse  Matrix  Reordering  Routines 


C***  Subroutine  ODRV 

C***  Driver  for  sparse  matrix  reordering  routines 
C 

SUBROUTINE  ODRV 

* (N,  U.JA.A,  P,IP,  NSP.ISP,  PATH,  FLAG) 


C 

C PA 

C Cl 

C 

C 

C 

C 

C 

c 

c 

C Class 

C 

C 

C 

C in 

C in^ 

C CO 

C in 

C ex 

C 
C 
C 
C 
C 

C wo 

C 

C 

C 

C 

C 

C 

C lA 

C lA 

C de 

C 

C tr 

C th 

C 

C vn 
C vfa 
C 
C 
C 

C vna 
C 

C vna 


PARAMETERS 

Class  abbreviations: 

V - supplies  a VALUE  to  the  driver 

r - contains  a RESULT  returned  by  the  driver 

1 - is  used  INTERNALly  by  the  driver 

a - is  an  ARRAY 

n - Is  an  INTEGER  variable 

f - Is  a REAL  variable. 


Parameter 


The  nonzero  entries  of  the  matrix  M are  stored  row  by  row 
In  the  array  A.  The  array  JA  contains  the  corresponding  column 
Indices;  l.e.  If  A(K)  - then  JA(K)  » J.  The  array  lA 

contains  pointers  to  delimit  the  rows  of  M;  IA(1)  Is  the  index 
in  JA  and  A of  the  first  entry  stored  In  the  Ith  row  of  M.  For 
example,  the  symmetric  5 by  S matrix 


would  be  stored  as 
112  3 4 


JA|1  34234455 
AI123456789 

IA(I+1)  - IA(I)  Is  the  number  of  nonzero  entries  in  row  I,  so 
IA(N+1),  where  N Is  the  number  of  rows  In  M,  Is  needed  to 
determine  the  length  of  the  Nth  row  In  A. 

If  M Is  a symmetric  matrix  only  the  elements  of  the  upper 
triangle  (Including  the  diagonal)  need  be  stored  In  A,  although 
the  full  matrix  may  be  stored. 


the  number  of  rows/columns  In  matrix  M. 
coefficient  matrix  for  the  system  of  linear  equations 
Mx  =•  b,  stored  in  compressed  form. 

Size  “ number  of  nonzeros  in  M (or  only  In  upper 
triangle  of  M for  symmetric  matrix), 
pointers  to  first  elements  of  each  row  in  A. 

Size  - N+l. 

the  column  numbers  corresponding  to  elements  of  A. 
Size  - size  of  A. 


r 


1 


c 

c 

c 

c 

c 

C vn 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 


rn 


A mlnlmuffl  degree  ordering  of  the  matrix  is  done  (ORDER).  If 
M is  symmetric  and  only  the  upper  triangle  is  being  stored,  the 
array  A needs  to  be  reordered. 


PATH  - information  on  %fhich  subroutines  are  to  be  called. 
Values  and  meanings  of  PATH  are: 

1 perform  minimum  degree  ordering  only. 

2 perform  minimum  degree  ordering  and 

reordering  of  symmetric  matrix.  This 
value  should  be  passed  only  if  M is 
symmetric  and  only  the  nonzeros  of  the 
upper  triangle  of  M are  being 
stored.  If  the  nonzeros  of  the  entire 
matrix  are  being  stored,  PATH  should 
equal  1. 

FLAG  - flag  for  error  return  from  subroutines.  Error  values 
and  their  meanings  are: 


0 

N+I 

9N+I 

lON+1 

llN+1 


no  error 

null  row  in  A — I 
ORDER  storage  exceeded  on  row  I 
ISP  too  small  to  allocate  space 
PATH  out  of  bounds 


c 

rna 

1 IP 

c 

1 

c 

rna 

I p 

fa 


The  result  of  the  ordering  algorithm  is  a permutation  of  the 
row  numbers  of  M and  the  inverse  of  the  permutation.  The  order 
of  the  columns  is  the  same  as  for  the  rows. 


- inverse  of  the  ordering  of  the  rows/columns  of  M. 

Size  = N. 

- ordering  of  the  rows/columns  of  M. 

Size  = N. 


Workspace  Is  needed  to  hold  the  temporary  vectors  used  in  the 
ordering  routine. 


NSP  - dimensioned  size  of  ISP.  NSP  must  generally  be  at 
least  2*(number  of  pairs  (I,J)  such  that  M(I,J) 
or  M(J,I)  is  nonzero)  + N for  the  ordering  routine 
and  at  least  N + size  of  A for  the  symmetric 
reordering  routine. 

ISP  - storage  space  divided  up  for  various  arrays  of 
the  subroutines. 


INTEGER  IA(1),  JA(1),  P(l),  IP(l), 

REAL  A(l),  ISP(l) 


PATH,  FLAG,  VV,  TMP,  Q 


IF  (PAXH.LT.  I .OR.  PATH.GT.  2)  GO  TO  HI 
C******  Allocate  space  for  ordering  subroutine  ********************** 
MAX  » NSP/2 
VV  = 1 

LV  - VV  + MAX 
IF  (MAX.LT.N)  GO  TO  110 

C******  Call  minimum  degree  ordering  routine  ************************ 
FLAG  - 0 
CALL  ORDER 

* (N,  lA,  JA,  P,  IP,  MAX,  ISP(VV),  ISP(LV),  FLAG) 

IF  (FLAG.NE.O)  GO  TO  100 


C******  Allocate  space,  call  symaetrlc  reorder  routine 
IF  (PATH.LT. 2)  CO  TO  1 
TMP  - 1 
Q - TMP  + N 

IF  (NSP+l-Q  .LT.  IA(N+1)-1)  GO  TO  110 
CALL  SRO 

* (N,  IP,  lA,  JA,  A,  ISP(TMP),  ISP(Q)) 

1 RETURN 

C 

C **  ERROR:  Error  Detected  In  ORDER 
100  RETURN 

C **  ERROR:  Insufficient  Storage 

110  FLAG  - 10*N  + I 

RETURN 

C **  ERROR:  Illegal  PATH  Specified 

111  FLAG  - 11*N  + 1 

RETURN 

END 

C 


C***  Subroutine  ORDER 

C**s  Minimum  degree  or^ering  algorithm  with  threshhold  search 
C 

SUBROUTINE  ORDER 

* (N,  IA,JA,  P,  IP,  MAX,  VV,LV,  FUG) 


C 

C 

C 

C 

C 

C nla 
C 
C 
C 

C nla 
C 

C nv 
C 


Input  variables:  N,  IA,JA,  MAX 

Output  variables:  P, IP,  FLAG 

Parameters  used  internally: 

I VV  - value  field  of  a linked  list  describing  adjacencies  of 
I vertices. 

I Size  .ge.  number  of  pairs  (I,J)  such  that  M(I,J)  or 

I M(J,I)  is  nonzero. 

1 LV  - link  field  of  the  linked  list. 

I Size  “ size  of  VV. 

I MAX  - dimensioned  size  of  VV  and  LV. 

INTEGER  U(l),  JA(l),  P(l),  IP(1),  VV(1),  LV(I),  FLAG, 

DMIN,  DTHR,  S,  SFS,  TMP,  VI,  VJ,  VK,  VL 


C * Initialize  free  storage  ****************************************** 
VK  - 1 

IF  (MAX.LT.N)  CO  TO  109 
DO  1 S-N,MAX 

1 LV(S)  - S+l 

LV(MAX)  - 0 

SFS  - 1 

C * Initialize  ordering,  degree,  and  adjacency  ************************ 
DO  2 K-1,N 
VK  - P(K) 

IP(VK)  - K 
VV(K)  - N+1 

2 LV(K)  - K 

SFS  - SFS  + N 


C * Initialize  nonzero  structure  ************************************** 
C * For  every  vertex  VK  *********************************************** 
C ***  For  every  vertex  VJ  adjacent  to  VK  ****************************** 
DO  8 VK-1,N 
JMIN  - lA(VK) 

JMAX  - IA(VK+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  101 
DO  7 J-JMIN,JMAX 
VJ  - JA(J) 

IF  (VJ.EQ.VK)  GO  TO  7 

C *****  Search  for  VJ  In  adjacency  of  VK  ****************************** 

LLK  - VK 

3 LK  - LLK 
LLK  - LV(LK) 

IF  (VV(LLK)  - VJ)  3,  7,  4 

C *****  Insert  VJ  In  adjacency  of  VK  ********************************** 

4 VV(VK)  = VV(VK)  + 1 

IF  (SFS.EQ.O)  GO  TO  109 
LLK  - SFS 
SFS  - LV(SFS) 

W(LLK)  =■  VJ 

LV(LLK)  - LV(LK) 

LV(LK)  - LLK 

C *****  Search  for  VK  In  adjacency  of  VJ  ****************************** 

LLJ  - VJ 

5 LJ  - LLJ 
LLJ  » LV(LJ) 

IF  (VVdlJ)  - VK)  5,  7,  6 

C *****  Insert  VK  In  adjacency  of  VJ  ********************************** 

6 VV(VJ)  « VV(VJ)  + 1 

IF  (SFS.EQ.O)  GO  TO  109 
LLJ  - SFS 
SFS  - LV(SFS) 

VV(LLJ)  = VK 

LV(LLJ)  - LV(LJ) 

LV(LJ)  =•  LLJ 

7 CONTINUE 

8 CONTINUE 
C 

C * Minimum  degree  algorithm  with  threshhold  search  ******************* 
C * Initialize  vertex  count  and  threshholds  for  search  **************** 
I - 0 
JMIN  - 1 
DTHR  - 0 
DMIN  - N+N 

C * While  unellmlnated  vertices  exist  ********************************* 
C ***  Search  for  vertex  VI  of  minimum  degree  ************************** 

9 JMIN  - MAXO  (JMIN,  I+l ) 

DO  10  J=JMIN,N 

VI  - P(J) 

IF  (VV(VI).LE.DTHR)  GO  TO  11 
10  DMIN  - MINO  (DMIN,  VV(VI)) 

JMIN  - 1 
DTHR  - DMIN 
DMIN  - N+N 
GO  TO  9 


C ***  Number  vertex  VI  of  minimum  degree  ****************************** 

11  JMIN  - J 
I - I+l 
VJ  - P(I) 

P(J)  - VJ 
IP(VJ)  - J 
P(I)  - VI 
IP(VI)  - I 
NI  - 1 

C ***  Delete  eliminated  vertices  from  adjacency  of  VI  ***************** 
C ***  For  every  vertex  VK  adjacent  to  VI  ****************************** 

LLI  - VI 

KHAX  - (VV(Vl)  - NI)  - N 
IF  (KHAX.LF..O)  00  TO  14 
DO  13  K-1,KMAX 

12  LI  - LLI 
LLI  - LV(LI) 

VK  - VV(LLI) 

C *****  If  VK  eliminated,  then  delete  from  adjacency  of  VI  ************ 
IF  (IP(VK).CT.I)  CO  TO  13 
LV(Ll)  - LV(LLI) 

LV(LLl)  - SFS 
SFS  - LLI 
LLI  - LI 
GO  TO  12 

13  CONTINUE 

C ***  Eliminate  vertex  VI  ********************************************* 
C ***  For  every  vertex  VK  ddj^cent  to  VI  ****************************** 

14  LLI  - VI 

KMAX  - (VV(VI)  - NI)  - N 
IF  (KMAX.LE.O)  GO  TO  21 
DO  20  K-1,KMAX 
LI  - LLI 
LLI  - LV(LI) 

VK  - VV(LLI) 

C *****  Merge  adjacency  of  VI  into  adjacency  of  VK  ******************** 
C *****  For  every  vertex  VJ  edjecent  to  VI  **************************** 
LLK  - VK 
LJ  - VI 

JMAX  =•  (VV(VI)  - NI)  - N 
IF  (JMAX.LE.O)  GO  TO  19 
DO  18  J=-1,JMAX 
LJ  = LV(LJ) 

VJ  =■  VV(LJ) 

IF  (VJ.EQ.VK)  GO  TO  18 

C *******  Search  for  VJ  in  adjacency  of  VK  ...  ************************ 

15  LK  - LLK 
LLK  - LV(LK) 

VL  = VV(LLK) 

IF  (VJ.LE.VL)  GO  TO  17 

C *******  ...  while  deleting  eliminated  vertices  ********************** 
IF  (IP(VL).GT.I)  GO  TO  16 
LVfLK)  - LV(LLK) 

LV(LLK)  - SFS 
SFS  - LLK 
LLK  - LK 
GO  TO  15 


16 


I 
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C *******  Insert  VJ  In  adjacency  of  VK  ******************************** 

17  IF  (VJ.EQ.VL)  GO  TO  18 

VV(VK)  - VV(VK)  + 1 
IF  (SFS.EQ.O)  CO  TO  109 
LLK  - SFS 
SFS  - LV(SFS) 

VVatK)  - VJ 
LV(LLK)  - LV(LK) 

LV(LK)  - LLK 

18  CONTINUE 

C *****  If  VK  of  minimal  degree,  then  number  vertex  VK,  ...  *********** 

19  IF  {VV(VK).GT.VV(V1))  GO  TO  20 

I - I+l 
J - IP(VK) 

VJ  - P(I) 

P(J)  - VJ 
IP(VJ)  - J 
P(I)  - VK 
IP(VK)  - I 
NI  - NI  + 1 

C *****  ...  recover  storage  for  adjacency  of  VK,  ********************** 

IMP  - LV(VK) 

LV(VK)  - SFS 
SFS  - TMP 

C *****  ...  and  delete  VK  from  adjacency  of  VI  ************************ 

LV(LI)  - LV(LLI) 

LV(LLI)  - SFS 
SFS  - LLI 
LLI  - LI 

20  CONTINUE 

C ***  Update  degrees  of  unellmlnated  vertices  adjacent  to  VI  ********** 
C ***  pof  every  vertex  VK  In  adjacency  of  VI  ************************** 

21  LI  - VI 

KMAX  •=  (VV(VI)  - NI)  - N 
IF  (KMAX.LE.O)  GO  TO  24 
DO  23  K-l,KMAX 
LI  - LV(LI) 

VK  - VV(LI) 

C *****  Update  degree  of  VK  and  threshholds  for  cyclic  search  ********* 
VV(VK)  - VV(VK)  - NI 
IF  (VV(VK).CE.DMIN)  CO  TO  23 
IF  (VV(VK).GT.DTHR)  GO  TO  22 
DMIN  - DTHR 
DTHR  - VV(VK) 

JMIN  = MINO  (JMIN,  IP(VK)) 

GO  TO  23 

22  DMIN  = VV(VK) 

23  CONTINUE 

C ***  Recover  storage  for  adjacency  of  VI  **************************** 

24  TMP  - LV(VI) 

LV(VI)  - SFS 
SFS  = TMP 

IF  (I.LT.N)  GO  TO  9 


C 


FLAG  - 0 
RETURN 


no  noonnooononnnonoonnnooo  onnonn  o on 


**  ERROR:  Null  row  In  A 
101  FLAG  - N + VK 
RETURN 

**  ERROR:  Insufficient  Storage 
109  FLAG  - 9*N  + VK 
RETURN 
END 


Subroutine  SRO 

Symmetric  reordering  of  sparse  symmetric  matrix 
SUBROUTINE  SRO 

* (N,  IP,  IA,JA,A,  TMP,  Q) 

Input  variables:  N,  IP,  IA,JA,A 

Output  variables:  IA,JA,A 


nla 


nla 


Parameters  used  Internally: 

I TMP  - Initially,  TMP(K)  Is  set  to  the  ntmber  of  elements 
I which  will  appear  In  the  Kth  row  of  M after 

I reordering.  Then  TMP  Is  Initialized  to  lA  and 

I USED  to  set  Q. 

I Size  = N. 

I Q - Initially,  Q(J)  Is  set  to  the  row  In  which  A(J) 

I (the  element  of  the  old  A)  will  appear  after 

I reordering.  Then  It  is  set  to  the  Index  of  A(J)  in 

I the  reordered  matrix. 

I Size  “ number  of  nonzeros  In  the  upper  trlaiijjle  of  M. 


The  subroutine  does  not  rearrange  the  order  of  the  rows,  but 
arranges  each  row  so  that  the  elements  which  will  be  above  the 
diagonal  after  reordering  are  filled  In.  If  M(1,J)  Is  above  the 
diagonal  but  Is  below  the  diagonal  after  reordering,  then  M(J,I) 
must  be  filled  In,  so  some  elements  will  appear  on  different  rows 
after  SRO  Is  finished. 


INTEGER  IP(1),  IA(l),  JA(1),  TMP(l),  Q(l),  QK 
REAL  A(l) 


******  Initialize  TMP  ******************************************** 
DO  1 I-1,N 

1 TMP(I)  =•  0 

Q AAAAA*  For  edch  row  of  A *4t**AAAA(AAA*A*AAA*4r*AAA*AA**AA*A****A**** 

DO  3 1=1,  N 
JMIN  = IA(I) 

JMAX  - IA(I+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  3 
Q AAAAAA  For  each  element  of  the  row 
DO  2 J-JMIN,JMAX 
K = JA(J) 

C ******  Adjust  TMP,  Q,  and  adjust  JA  if  necessary  ***************** 
IF  (IPfKl.LT.Ii^d))  JA(J)  - I 
IF  (IP(K).GE.IP(I))  K - I 
Q(J)  » K 

2 TMP(K)  - TMP(K)  + 1 

3 CONTINUE 


o o 


C ******  Set  new  lA  and  copy  It  Into  TMP  *************************** 
DO  4 I-1,N 

U(I+l)  - IA(I)  + TMP(I) 

4 TMP(I)  - IA(I) 

C ******  Q(j)  gets  position  of  A(J)  after  reordering  *************** 
JMIN  - IA(1) 

JMAX  - 1A(N+1)  - 1 
DO  5 J-JMIN,JMAX 
K - Q(J) 

Q(J)  - TMP(K) 

5 TMP(K)  - TMP(K)  + 1 


******  Reset  JA  and  A ******************************************** 

DO  7 J=JMIN,JMAX 

6 IF  (Q(J).Eg.J)  GO  TO  7 

K - Q(J) 

Q(J)  - Q(K) 

Q(K)  - K 
JAK  - JA(K) 

JA(K)  - JA(J) 

JA(J)  - JAK 
AK  - A(K) 

A(K)  - A(J) 

A(J)  - AK 
GO  TO  6 

7 CONTINUE 
RETURN 
END 
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Subroutines  for  Solving  Sparse  Symmetric  PoBltlve 
Definite  Systems  of  Linear  Equations 


C***  Subroutine  SDRV 

C***  Driver  for  subroutines  for  solving  sparse  symmetric  positive 
C definite  systems  of  linear  equations 

C 

SUBROUTINE  SDRV 

* (N,  P,IP,  lA.JA.A,  B,  Z,  NSP,ISP,RSP,  PATH,  FLAG) 

C 

C PARAMETERS 

C Class  abbreviations  ate: 

C v - supplies  a VALUE  to  the  driver 

C r - contains  a RESULT  returned  by  the  driver 

C 1 - is  used  INTERNALly  by  the  driver 

C a - is  an  ARRAY 

C n - Is  an  INTEGER  variable 

C f - Is  a REAL  variable. 

C 

C Class  I Parameter 
C I 

C The  nonzero  entries  of  the  matrix  M are  stored  row  by  row 

C in  the  array  A.  The  array  JA  contains  the  corresponding  coition 

C indices;  l.e.  if  A(K)  - M(1,J),  then  JA(K)  - J.  The  array  U 

C contains  pointers  to  delimit  the  rows  of  M — IA(1)  Is  the  Index 

C in  JA  and  A of  the  first  entry  stored  In  the  Ith  row  of  M. 

C Only  the  nonzero  entries  on  or  above  the  diagonal  need  be  stored. 

C However,  the  subroutines  will  work  if  all  nonzeros  are  stored. 

C For  example,  the  symmetric  5 by  5 matrix 

C 1 0 2 3 0 

C 0 4 0 0 0 

C 2 0 5 6 0 

C 3 0 6 7 8 

C 0 0 0 8 9 

C would  be  stored  as 

C 1123456789 


IA(I+l)  - IA(I)  Is  the  number  of  nonzero  entries  In  row  I,  so 
IA(N+1),  where  N is  the  number  of  rows  In  M,  is  needed  to 
determine  the  length  of  the  Nth  row  In  A. 


C vfa 
C 
C 
C 

C vna 
C 

C vna 
C 

C vna 


- the  number  of  rows/columns  In  matrix  M. 

- coefficient  matrix  for  the  system  of  linear  equations 

Mx  = b,  stored  in  compressed  form. 

Size  " number  of  nonzeros  in  upper  triangle  of  M 
(or  the  number  of  nonzeros  In  all  of  M) . 

- pointers  to  the  first  element  of  each  row  In  A. 

Size  = N+1. 

- the  column  numbers  corresponding  to  elements  of  A. 

Size  “ size  of  A. 

- right-hand  side  for  the  equation  Mx  » b.  B and  Z 

cannot  be  the  same  vector. 

Size  - N. 


o n 


C rna 
C 


solution  vector  for  the  equation 
cannot  be  the  same  vector. 

Size  “ N. 


R and  1. 


The  solution  of  the  system  is  done  in  three  stages: 

SYMFAC  - The  matrix  M is  processed  symbolically  to  determine 
where  fillin  will  occur  during  factorization. 

NUMFAC  - The  matrix  M is  factored  numerically  into  two 
triangular  matrices. 

NUMSLV  - The  system  resulting  from  NUMFAC  is  solved. 

For  several  systems  with  identical  nonzero  structures,  SYMFAC 
need  be  done  only  once,  then  NUMFAC  and  NUMSLV  are  done  for  each 
system.  For  several  system  with  identical  matrices  M and 
different  right-hand  sides,  SYMFAC  and  NUMFAC  need  be  done  only 
once,  then  NUMSLV  is  done  for  each  right-hand  side. 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C vna 


PATH  - information  on  which  subroutines  are  to  be  called. 
Values  and  meanings  of  PATH  arc: 

1 perform  SYMFAC,  NUMFAC  and  NUMSLV. 

2 perform  NUMFAC  and  NUMSLV.  (SYMFAC  is 

assumed  to  have  been  done  in  a manner 
compatible  with  the  driver's  storage 
allocation.) 

3 perform  NUMSLV  only.  (SYMFAC  and  NUMFAC 

are  assumed  to  have  been  done.) 

FLAG  - flag  for  error  return  from  subroutines.  Error  values 
and  their  meanings  are: 

0 no  error 

N+I  row  I of  A is  null 

2N+I  duplicate  entry  on  row  1 of  A 
6N+I  storage  exceeded  on  row  I in  SYMFAC 
TN+l  storage  exceeded  in  NUMFAC 

8N+I  diagonal  elemenL“0  on  cow  1 In  NUMFAC 

lON+1  ISP/RSP  too  small  to  allocate  space 

llN+l  PATH  out  of  bounds 


The  rows  and  columns  of  the  original  matrix 
arbitrarily  reordered  before  calling  the  driver. 


M can  be 

If  no  reordering 


is  done,  then  P(I)  = IP(I)  » I for 
is  returned  in  the  original  order. 


The  answer  vector 


the  ordering  of  the  rows  (and  columns)  of  M.  P(I) 
is  the  number  of  the  row  of  M which  becomes  the 
Tth  row  after  reordering. 

Size  = N. 

the  inverse  of  the  ordering  of  the  rows  of  M.  That 
is,  IP(P(I))  = I for  I-l,N. 

Size  = N. 


Workspace  is  needed  to  hold  the  factored  form  of  the  matrix 
plus  various  temporary  vectors. 

1 ISP  - integer  storage  space  divided  up  for  various  arrays 
I of  the  subroutines.  ISP  and  RSP  should  be  the 

I same  array.  This  allows  declaration  of  all  real 

I storage  to  be  double  precision. 

I NSP  - dimensioned  size  of  ISP  and  RSP.  NSP  generally 
I must  be  at  least  4N+1  + 2*K  (where  K « (number  of 

1 nonzeros  in  the  upper  triangle  of  M) ) , since  ISP 


f 


c 

c 

c 

c 

C fa 

C 

C 

C 

C 
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c 

c 

c****** 
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I and  RSP  must  hold: 

I four  vectors  of  fixed  length; 

I JU  (with  size  - K + flllln  - compression); 

I IJ  (with  size  » K + flllln). 

I RSP  - real  storage  space  divided  up  for  various  arrays  of 
I the  subroutines.  ISP  and  RSP  should  be  the  same 

I array.  This  allows  declaration  of  all  real  storage 

I to  be  double  precision. 

INTEGER  P(l),  IP(1),  IA(1),  JA(1),  ISP(l),  PATH,  FLAG, 

Q,  D,  U,  ROW.  TMP,  UMAX 
REAL  A(l),  B(l),  Z(l),  RSP(l) 

EQUIVALENCE  (ISP(l),  RSP(l)) 

IF  (PATH.LT.l  .OR.  PATH.GT.3)  GO  TO  111 
Initialize  and  divide  up  temporary  st<)rage  ****************** 


IJU 

- 

1 

lU 

IJU 

+ 

N 

JL 

■ 

lU 

+ 

N+1 

JU 

- 

JL 

+ 

N 

Q 

= 

NSP 

- 

N 

JUMAX  - Q - JU 

IF  (JUMAX. LT.O)  CO  TO  110 


Q*Aik**A  Cdll  Sllbroutin6S  ****A****4t******A*****AA***A***A*4t***AAA***4r 

FLAG  - 0 

IF  (PATH.CT. 1)  GO  TO  1 
CALL  SSF 

* (N,  P.  IP,  lA,  JA,  ISP(IJU),  ISP(JU),  ISP(IU),  JUMAX, 

* RSP(Q),  ISP(JL),  FLAG) 

IF(FLAG.NE.O)  GO  TO  100 

C 

1 n = Q - N 

U = JU  + ISP(IJU+(N-1 )) 

ROW  » Q 
UMAX  - D - U 
IF  (PATH.CT. 2)  CO  TO  2 
CALL  SNF 

* (N,  P.  IP,  lA,  JA,  A, 

* RSP(D),  ISP(IJU),  ISP(JU),  ISP(IU),  RSP(U),  UMAX, 

* RSP(ROW),  ISP(JL),  FLAG) 

IF  (FLAG.NE.O)  GO  TO  100 

C 

2 TMP  - Q 
CALL  SNS 

* (N,  P,  RSP(D),  ISP(IJU),  ISP(JU),  ISP(IU),  RSP(U),  Z,  B, 

* RSP(TMP)) 

RETURN 

C 

C **  ERROR:  Error  Detected  in  SSF,  SNF,  or  SNS 
100  RETURN 

C **  ERROR:  Insufficient  Storage 
no  FLAG  = 10*N  + 1 

RETURN 

C **  ERROR:  Illegal  PATH  Specification 
111  FLAG  = 11*N  + 1 

RETURN 
END 
C 

C 


n n n 
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C YALE  SPARSE  MATRIX  PACKAGE  - SYMMETRIC  CODES 

C SOLVING  THE  SYSTEM  OF  EQUATIONS  Mx  - b 

C 

C I.  SUBROUTINE  NAMES 

C Subroutine  names  are  of  the  form  Sxx  where: 

C (1)  the  first  letter  Is  S for  symmetric  matrices; 

C (2)  the  second  letter  Is  either  S for  symbolic  or  N for 

C numerical  processing; 

C (3)  the  third  letter  Is  either  F for  factorization  or  S for 
C solution. 

C 

C II.  CALLING  SEQUENCES 

C The  Input  matrix  can  be  processed  with  an  ordering  subroutine 

C before  using  the  remaining  subroutines.  If  this  Is  done  and  only 
C the  upper  triangle  of  M Is  being  stored,  SRO  should  be  called  to 
C reorder  the  matrix  Into  symmetric  form  before  using  the  other 
C subroutines.  If  an  ordering  subroutine  Is  not  used,  set  P(I)  ~ 

C IP(I)  = I for  1=1, N.  Then  the  calling  sequence  Is 
C SSF  (symbolic  factorization) 

C SNF  (numerical  factorization) 

C SNS  (called  once  for  each  right-hand  side). 

C 

C III.  STORAGE  OF  SPARSE  MATRICES 

C The  nonzero  entries  of  the  matrix  M are  stored  row  by  row 

C In  the  array  A.  The  array  JA  contains  the  corresponding  column 

C Indices;  l.e.  If  A(K)  - M(I,J),  then  JA(K)  - J.  The  array  lA 

C contains  pointers  to  delimit  the  rows  of  M — IA(I)  Is  the  Index 

C In  JA  and  A of  the  first  entry  stored  In  the  Ith  row  of  M. 

C Only  the  nonzero  entries  on  or  above  the  diagonal  need  be  stored. 

C However,  the  subroutines  will  work  If  all  nonzeros  are  stored. 

C For  example,  the  symmetric  5 by  5 matrix 
C 1 0 2 3 0 

C 0 4 0 0 0 

C 2 0 5 6 0 

C 3 0 6 7 8 

C 0 0 0 8 9 

C would  be  stored  as 

C 1123456789 

C lA  I 1 4 5 7 9 10 

C JA  1134234455 

C AI123456789 

C 

C IA(I+1)  - IA(I)  Is  the  number  of  nonzero  entries  In  row  I,  so 

C IA(N+1),  where  N Is  the  number  of  rows  In  M,  Is  needed  to 

C determine  the  length  of  the  Nth  row  In  A. 

C The  unit  triangular  matrix  U Is  stored  In  a similar  fashion 

C using  the  arrays  lU,  JU,  and  U except  that  an  additional  vector 
C IJU  Is  used  to  compress  storage  of  JU.  IJU(K)  points  to  the 

C starting  location  In  .lU  of  entries  for  the  Kth  row.  Compression 

C occurs  In  two  ways.  First,  If  a row  1 was  merged  Into  the  current 
C row  K,  and  the  number  of  elements  merged  In  from  row  I (some  tall 

C portion  of  row  I)  Is  the  same  as  the  final  length  of  row  K,  then 

C the  Kth  row  and  the  tall  are  Identical  and  IJU(K)  can  point  to 

the  start  of  the  tall.  Second,  If  some  tall  portion  of  the  K-lst 
row  equals  the  head  of  the  Kth  row,  then  IJU(K)  can  point  to  the 
start  of  that  tall  section.  For  example,  the  nonzero  structure  of 
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C the  matrix 

C d 0 X X X 

C 0 d 0 X X 

C 0 0 d X 0 

C 0 0 0 d X 

C 0 0 0 0 d 

C might  be  stored.  Ignoring  the  diagonal,  as 
C 11  2 3 4 5 6 

C lU  I 1 4 6 7 8 8 

C JU  I 3 4 5 4 

C IJU  I 1 2 4 3 

C 

C The  diagonal  entries  of  U are  assumed  to  be  equal  to  one 

C and  are  not  stored.  The  array  D contains  reciprocals  of  entries 
C of  the  diagonal  matrix  In  the  U DU  decomposition. 

C 

C IV.  ADDITIONAL  STORAGE  SAVINGS 

C In  SSF  and  SNF,  P and  IP  can  be  the  same  vector  In  the 

C calling  sequences  If  no  reordering  of  the  matrix  has  been  done 

(l.e.  P(I)  - IP(I)  - 1 for  I-1,N). 

In  SNS,  B ard  Z can  be  the  same;  however,  the  right-hand 
side  B will  be  destroyed. 

V.  PARAMETEHS 

Following  Is  a list  of  parameters  to  the  programs.  Names  are 
uniform  among  the  various  subroutines.  Class  abbreviations  are: 

V - supplies  a VALUE  to  the  subroutine 
r - contains  a RESULT  returned  by  the  subroutine 
1 - Is  used  INTERNALly  by  the  subroutine 
a - Is  an  ARRAY 
n - Is  an  INTEGER  variable 
f - Is  a REAL  variable. 


Class  I Parameter 


fva 


fva 

fvra 

nvra 

nr 


C 

C 

C 

C 

C nvra 

C 

C 

C nva 
C 
C 
C 


A - coefficient  matrix  for  the  system  of  linear  equations 
^bc  > b,  stored  In  compressed  form. 

Size  “ either  the  number  of  nonzeros  In  the  upper 

triangle  of  M,  or  the  number  of  nonzeros  In 
all  of  M (see  section  III). 

B - right-hand  side  for  the  equation  Mx  ■ b. 

Size  - N. 

D - Inverse  of  diagonal  matrix  In  UtDU  factorization 
(also  used  for  temporary  results  In  SNF). 

Size  - N. 

lA  - pointers  to  first  elements  of  each  row  In  A. 

Size  - N+1. 

FLAG  - flag  for  error  return  from  subroutines.  Error  values 
and  their  meanings  are: 

0 no  error 

N+I  row  I of  A Is  null 

2N+I  duplicate  entry  on  row  I of  A 
6N-M  JU  storage  exceeded  on  row  I 
7N+1  U storage  exceeded 

8N-M  zero  diagonal  element  on  row  I 

IJU  - pointers  to  the  first  elements  of  each  row  In  JU, 
used  to  compress  storage  of  JU. 

Size  - N. 

IP  - Inverse  of  the  ordering  of  the  rows  of  M.  For 

example.  If  row  1 Is  the  5th  row  after  reordering, 
then  IP(l)-5. 

Size  « N. 


i 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


nvra 

nvra 

nvra 

nv 

nv 

nva 


fvra 


nv 

fra 


lU  - pointers  to  the  first  elements  of  each  row  in  U. 
Size  « N-fl. 

JA  - column  numbers  corresponding  to  elements  of  A. 

Size  - size  of  A. 

JU  - column  numbers  corresponding  to  elements  of  U. 

Size  - size  of  U - compression. 

JUMAX  - declared  dimension  of  JU. 

N - nunber  of  rows/columns  in  matrix  M. 

P - ordering  of  rows  (and  columns)  of  M.  P(I)  is 

the  number  of  the  row  of  M which  becomes  the  1th 
row  after  reordering. 

Size  - N. 

U - upper  triangular  matrix  resulting  from  the 

factorization  of  M,  stored  in  compressed  form. 
Size  - number  of  nonzeros  in  upper  triangle  of  M 
plus  flllln  (IU(N+1)-1  after  SSF) . 

UMAX  - declared  dimension  of  U. 

Z - solution  vector  for  the  equation  Mx  « b. 

Size  - N. 


C***  Subroutine  SSF 

C***  Symbolic  Ut-I)-U  factorization  of  sparse  symmetric  matrix 
C 

SUBROUTINE  SSF 

* (N,  P,IP,  1A,JA,  1JU,JU,IU, JUMAX,  Q,  JL,  FLAG) 


nla 


nla 


Input  variables: 
Output  variables: 


N,  P,IP,  lA.JA,  JUMAX 
1JU,JU,IU,  FLAG 


Parameters  used  Internally: 

JL  - linked  list  of  rows  to  be  merged.  IF  the  Kth  row  is 
being  processed,  JL(K)  contains  the  number  of  the 
first  row  to  be  merged  with  the  Kth  row,  JL(JL(K)) 
is  the  number  of  the  second  row,  etc. 

Size  “ N. 

Q - Suppose  M'  is  the  result  of  reordering  M.  If 

processing  of  the  Kth  row  of  M'  (hence  the  Kth  row 
of  U)  is  being  done,  Q(J)  is  initially  nonzero  if 
M'(K,J)  is  nonzero  and  above  the  diagonal.  Since 
values  need  not  be  stored,  each  entry  points  to  the 
next  nonzero  and  Q(K)  points  to  the  first.  N+1 


For  example,  if  N-9  and 


indicates  the  last  element, 
the  5th  row  of  M'  is 
OxxOxOOxO 
then  Q will  initially  be 

aaaaSaalOa  (a-  arbitrary). 

As  the  algorithm  proceeds,  other  elements  of  Q 
Inserted  in  the  list  because  of  flllln. 

Size  “ N. 


Internal  variables: 

JUMIN,JUPTR  - are  the  indices  in  JU  of  the  first  and  last 
elements  in  either  the  last  or  the  current  row. 
LMAX  - length  of  longest  row  merged  into  Q. 

LUI  - number  of  elements  in  a row  to  be  merged  into  Q. 

LUK  - number  of  elements  in  the  current  row  (Q). 


are 


an  n n n o o 


INTEGER  P(l).  IP(l),  IA(l),  JA(l),  IJU(l),  JU(1),  IU(1), 

* Q(l),  JL(l),  FLAG,  VJ,  qM 

Initialize  ***A****ik*A******A******A*******A*************** 

JUMIN  - 1 
JUPTR  - 0 
IU(1)  - 1 
DO  I K-l,N 

1 JL(K)  - 0 

******  For  each  row  ********************************************** 

DO  15  K-1,N 

******  Initialize  Q to  structure  of  Kth  row  above  diagonal  ******* 
LUK  - 0 
0(K)  - N+1 
JMIN  - IA(P(K)) 

JMAX  - IA(P(K)+1)  - I 
IF  ( JMIN. GT. JMAX)  GO  TO  101 
DO  3 J-JMIN,JMAX 
VJ  - IP(JA(J)) 

IF  (VJ.LE.K)  GO  TO  3 
QM  - K 

2 M - QM 
QM  = Q(M) 

IF  (QM.LT.VJ)  GO  TO  2 
IF  (QM.EQ.VJ)  GO  TO  102 
LUK  - LUK+l 
Q(M)  • VJ 
Q(VJ)  «=  QM 

3 CONTINUE 


******  Compute  fillin  for  Q by  *********************************** 

LMAX  = 0 
IJU(K)  - JUPTR 
I - K 

Q Linking  through  JL  and  **********************ik************* 

A I - JL(I) 

IF  (I.EQ.O)  GO  TO  8 

LUI  » IU(I+1)  - (IU(I)+1) 

JMIN  = IJU(I)  + 1 

JMAX  = IJU(I)  + LUI 
IF  (LUI. LE. LMAX)  GO  TO  5 
LMAX  - LUI 
IJU(K)  - JMIN 

5 l)M  - K 

C ******  Merging  each  row  with  Q *********************************** 

DO  7 J-JMIN,JMAX 
VJ  - JU(J) 

6 M - QM 
QM  - Q(M) 

IF  (QM.LT.VJ)  GO  TO  6 
IF  (QM.EQ.VJ)  GO  TO  7 
LUK  - LUK+1 
Q(M)  - VJ 
Q(VJ)  - QM 
C^I  - VJ 

7 CONTINUE 
GO  TO  4 


****** 

8 

****** 


10 


11 


Check  If  row  duplicates  another 
IF  (LUK. EQ.LMAX)  CO  TO  14 
see  If  tall  of  K-lst  row  matches  head 
IF  (JUMIN.GT. JUPTR)  GO  TO  12 
I - Q(K) 

DO  9 JMIN-JUMIN, JUPTR 

IF  (JU(JMIN)-I)  9,  10,  12 
CONTINUE 
GO  TO  12 
IJU(K)  - JMIN 
DO  11  J-JMIN, JUPTR 

IF  (JU(J).NE.I)  GO  TO  12 
I - Q(I) 

IF  (I.GT.N)  GO  TO  14 
CONTINUE 
JUPTR  - JMIN  - 1 


If  not  ****************** 


of  Kth  ************** 


******  Set  Kth  row  of  U to  Q ************************************* 

12 


JUMIN  - JUPTR  + 1 

JUPTR  - JUPTR  + LUK 
IF  (JUPTR.GT. JUMAX) 
I - K 

DO  13  J -JUMIN, JUPTR 
I - Q(I) 

13  JU(J)  - I 

IJU(K)  - JUMIN 


GO  TO  106 


******  If  more  than  one  element  In  row,  adjust  JL 

14  IF  (LUK.LE.l)  CO  TO  15 

I - JU(IJU(K)) 

JL(K)  = JL(I) 

JL(I)  - K 

15  IU(K+1)  - IU(K)  + LUK 


FLAG  - 0 
RETURN 

**  ERROR:  Null  Row  In  A 


101  FLAG  - N + P(K) 

RETURN 

' **  ERROR:  Duplicate  Entry  In  A 

102  FLAG  - 2*N  + P(K) 

RETURN 

: **  ERROR:  Insufficient  Storage  for  JU 
106  FLAG  - 6*N  + K 
RETURN 
END 


C***  Subroutine  SNF 

C***  Numerical  Ut-D-U  factorization  of  sparse  symmetric  positive 


definite  matrix 


C 

C 

C 

C 

C 

C 

C nlva 

C 

C 

C 

C 

C 

C nla 

C 

C 

C nla 


SUBROUTINE  SNF 

(N,  P,IP.  lA.JA.A,  D.  lJU,  JU.IU.U.UMAX,  IL,  JL,  FLAG) 


Input  variables:  N,  P,IP,  1A,JA,A,  IJU,JU,IU 

Output  variables:  D,U,  FLAG 

Parameters  used  Internally: 

I  D - If  the  Kth  row  of  U Is  being  computed,  D(l)  through 
1 D(lC-l)  contain  reciprocals  of  the  entries  of  the 

I diagonal  matrix  D from  the  decomposition.  The 

I remainder  of  D Is  Initialized  to  the  structure  of 

I the  Kth  row  of  M (after  reordering)  and  Is  adjusted 

I to  become  the  Kth  row  of  U. 

I IL  - IL(I)  points  to  the  first  element  of  the  Ith  row  to  be 
I used  in  adjusting  the  current  row. 

I Size  = N. 

I JL  - linked  list  of  rows  to  be  used  tu  ailjustlng  the  current 
I row.  If  the  Kth  row  is  being  processed,  JL(K) 

I contains  the  number  of  the  first  row  to  be  used  with 

I the  Kth  row,  JL(JL(K))  is  the  number  of  the  second 

I row,  etc. 

I Size  = N. 


P(l),  IP(1),  IA(1),  JA(l),  IJU(l),  JU(1),  IU(1), 

IL(1),  JL(1),  FLAG,  VK,  VJ 


INTEGER  P(l),  IP(1),  IA(1) 

UMAX,  IL(1),  JL(1),  FL 

DIMENSION  A(l),  D(l),  U(l) 


******  Initialize  .IL,  check  storage  ** 
IF  (IU(N+1)-1  .GT.  UMAX)  GO  TO  107 
DO  1 K-1,N 
1 JI,(K)  = 0 


AAAAAA  for  each  row  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

DO  10  K=1,N 

******  Initialize  D on  and  above  the  diagonal  ******************** 
JMIN  = IU(K) 

UMAX  - IU(K+1)  - 1 
IF  (JMIN. GT. UMAX)  GO  TO  3 
MU  = IJU(K)  - IU(K) 

DO  2 J=JMIN,.1MAX 

2 D(JU(MU+J))  - 0 

3 D(K)  - 0 
VK  - P(K) 

JMIN  - lA(VK) 

JMAX  - IA(VK+1)  - 1 
DO  4 J-JMIN,JMAX 
VJ  - IP(JA(J)) 

IF  (K.LE.VJ)  D(VJ)  - A(J) 

4 CONTINUE 


For  each  element  In  lower  triangle  to  be  eliminated  ******* 
DK  - D(K) 

NXTI  - JL(K) 

I - NXTI 

IF  (I.EQ.O)  GO  TO  8 

Change  D and  adjust  IL  and  JL  ***************************** 
NXTI  - JL(I) 

UKIDI  - - U(IL(I))  * D(I) 

DK  = DK  + UKIDI  * U(IL{I)) 

U(1L(I))  - UKIDI 
JMIN  = IL(  I ) + 1 
JMAX  = IU(I+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  7 
MU  - IJU(I)  - IU(I) 

DO  6 J-JMIN,JMAX 

D(JU(MU+J))  - D(JU(MU+J))  + UKIDI  * U(J) 

IL(I)  - JMIN 
J - JU(MU+JMIN) 

JL(I)  = JL(J) 

JL{J)  = I 
GO  TO  5 


****** 

8 


10 


Set  D(K)  and  copy  rest  of  D Into  Kth  row  of  U 
IF  (DK.EQ.O)  GO  TO  108 
D(K)  = 1 / DK 
JMIN  - IU(K) 

JMAX  = IU(K+1)  - 1 
IF  (JMIN. GT. JMAX)  GO  TO  10 
MU  = IJU(K)  - JMIN 
DO  9 J=JMIN;JMAX 
U(J)  = D(JU(MU+J)) 

IL(K)  = JMIN 
I = JU(MU+-JMIN) 

JL(K)  = JL(I) 

JL(I)  = K 
CONTINUE 


************* 


FLAG  => 
RETURN 


C **  ERROR:  Insufficient  Storage  for 

107  FLAG  - 7*N  + 1 
RETURN 

C **  ERROR:  Zero  Pivot 

108  FLAG  - 8*N  + K 
RETURN 

END 


c 

c 

c 


C***  Subroutine  SNS 

C***  Numerical  solution  of  sparse  symmetric  positive  definite  system  of 
C linear  equations  given  Ut-D-U  factorization 

C 

SUBROUTINE  SNS 

* (N.  P,  D,  IJU.JU.IU.U,  Z,  B,  TMP) 

C 

C Input  variables:  N,  P,  D,  1JU,JU,U,  B 

C Output  variables:  Z 

C 

C Parameters  used  Internally: 

C fla  I TMP  - vector  which  gets  result  of  solving  Ut  Dy  « b. 

C I Size  - N. 

C 

INTEGER  P(l),  IJU(l),  JU(1),  IU(1) 

REAL  D(l).  U(l),  Z(l),  B(l),  TMP(l) 

C 

C ******  Initialize  TMP  to  the  reordered  B ************************* 
DO  1 K-l.N 

1 TMP(K)  = B(P(K)) 

C ******  Solve  Ut  Dy  * b by  forward  substitution  ***************** 
DO  3 K“1,N 
TMPK  = TMP(K) 

JMIN  - IU(K) 

JMAX  =■  lUCK+I)  - I 

IF  (JMIN. GT. JMAX)  GO  TO  3 
MU  = IJU(K)  - JMIN 
DO  2 J-JMIN,JMAX 

2 TMP(JU(MU+J))  - TMP(JU(MU+J))  + U (J ) * TMPK 

3 TMP(K)  = TMPK  * D(K) 

C 

C ******  Solve  Ux  = y by  back  substitution  *********************** 
K = N 

DO  6 1=1,  N 
SUM  = TMP(K) 

JMIN  = IU(K) 

JMAX  - IU(K+1)  - 1 
IF  (JMIN. GT,  JMAX)  GO  TO  5 
MU  = IJU(K)  - JMIN 
DO  4 J-JMIN,JMAX 

4 SUM  = SUM  + U(J)  * TMP(JU(MU+J)) 

5 TMP(K)  - SUM 
Z(P(K))  - SUM 

6 K - K-1 
RETURN 


n n 


I 

i 


Test  Driver  for  Sparse  Symmetric  Matrix  PackaKe 


C***  Program  STST 

C***  Test  Driver  for  Symmetric  Codes  In  Yale  Sparse  Matrix  Package 
C 

C Variables: 

C 

C NG  - size  of  grid  used  to  generate  test  problem. 

C 

C N - number  of  variables  and  equations  (■  NG  x NG). 

C 

C lA  - INTEGER  one-dimensional  array  used  to  store  row  pointers 

C to  JA  and  A;  DIMENSION  « N+1. 

C 

C JA  - INTEGER  one-dlmenslonal  array  used  to  store  column 

Indices  of  nonzero  elements  of  (upper  triangle  of)  M; 

DIMENSION  “ number  of  nonzero  entries  In  (upper  triangle 

C of)  M. 

C 

C A - REAL  one-dlmenslonal  array  used  to  store  nonzero  elements 

C of  (upper  triangle  of)  M;  DIMENSION  = number  of  nonzero 

C entries  in  (upper  triangle  of)  M. 

C 

C X - REAL  one-dlmenslonal  array  used  to  store  solution  x; 

C DIMENSION  = N. 

C 

C B - REIAL  one-dlmenslonal  array  used  to  store  right-hand-side  b 

C DIMENSION  = N. 

C 

C P - INTEGER  one-dlmenslonal  array  used  to  store  permutation  of 

C rows  and  columns  for  reordering  linear  system; 

C DIMENSION  = N. 

C 

C IP  - INTEGER  one-dimensional  array  used  to  store  Inverse  of 

C permutation  stored  In  P;  DIMENSION  - N. 

C 

C NSP  - declared  dimension  of  one-dlmenslonal  arrays  ISP  and  RSP. 

C 

C ISP  - INTEGER  one-dimensional  array  used  as  working  storage 

C (equivalenced  to  RSP);  DIMENSION  = NSP. 

C 

C RSP  - REAL  one-dlmenslonal  array  used  as  working  storage 

C (equivalenced  to  ISP);  DIMENSION  - NSP. 

C 

C 

INTEGER  lA(lOl),  JA(500),  P(IOO),  IP(IOO),  ISP(1500), 

* CASE,  PATH,  FLAG,  APTR,VP,VQ,  X,XMIN,XMAX,  Y,YMIN,YMAX 
REAL  A(500),  Z(IOO),  B(IOO),  RSP(1500) 

EQUIVALENCE  (ISP(l),  RSP(l)) 

DATA  NSP/1500/,  EPS/lE-5/ 

C 

INDEX  (I,  J)  » NG*I  + J - NG 
C 

NG  - 3 
N - NG*NG 


c****** 

c 


(^AAAAAA 


c 

QAAAAAA 


c 

(]  AAAAA 


c 

c****** 


For  CASE-1  ue  store  the  entire  matrix,  for  CASE-2  we  store 
only  the  upper  triangular  part 
DO  5 CASE-1,2 

Set  up  matrix  for  five  point  finite  difference  operator  ****** 
APTR  - 1 
DO  2 I-1,NG 
DO  2 J-1,NG 

VP  - INDEX  (I,  J) 

P(VP)  - VP 
IP(VP)  - VP 
lA(VP)  - APTR 
SUM  - 0 

XMIN  = MAXO  ( 1,  I-l) 

XMAX  - MINO  (NG,  I+l) 

YMIN  = MAXO  ( 1,  I-l) 

YMAX  = MINO  (NG,  J+1 ) 

DO  1 X^MIN,XMAX 
TX)  1 Y-YMIN,YMAX 

IF  ((X-I)  * (Y-J)  .NE.  0)  GO  TO  1 
VQ  - INDEX  (X,  Y) 

JA(APTR)  = VQ 
A(APTR)  - 4 

IF  (VP  .NE.  VQ)  A(APTR)  = -1 
SUM  = SUM  + A(APTR)  * VO 

If  CASE-2,  do  not  store  elements  below  diagonal  ************** 
IF(VP.GT.VQ  .AND,  CASE.EQ.2)  GO  TO  1 
APTR  = APTR  + 1 
CONTINUE 
B(VP)  = SUM 
CONTINUE 
IA(N+1)  = APTR 
NZA  » IA(N+1)  - 1 

Output  original  array  A ******** 

IF  (CASE.EQ. 1)  PRINT  1001,  NG,NG 

FORMAT  (/'  ***  FIVE-POINT  OPERATOR  ON  ',  II,  ' BY  ' II,  ' GRID  ' 
/ ' (ALL  ENTRIES  OF  MATRIX  STORED)  ') 

IF  (CASE.EQ.2)  PRINT  1002,  NG,NG 

FORMAT  (/'  ***  FIVE-POINT  OPERATOR  ON  ',11,  ' BY  ' II,  ' GRID  ' 
/ ' (ONLY  ENTRIES  OF  UPPER  nclANGLE  STORED)  ') 

PRINT  1003,  (IA(I),I=1,N),  IA(N+1 ) 

FORMAT  (/'  COEFFICIENT  MATRIX:  '/ 

/'  LA  (INDICES  OF  FIRST  ELEMENTS  IN  ROWS)' 

/(10I5)) 

PRINT  1004,  (I, JA(I),A(I),  1=1, NZA) 

FORMAT  (/'  JA  A ' 

/'  I COLUMN  INDICES  MATRIX' 

/(I3,  110,  F16.5)) 

PRINT  1005,  (B(I),  1=1, N) 

FORMAT  (/'  RIGHT  HAND  SIDE  B:  ' 

/(5F10.5)) 

Call  ODRV  *************************************************** * 
FLAG  = 0 
PATH  = CASE 
CALL  ODRV 

(N,  IA,JA,A,  P,IP,  NSP,RSP,  PATH,  FLAG) 

IF  (FLAG.NE.O)  GO  TO  101 

Output  reordered  array  A ************************************* 
PRINT  1006,  (I,P(I),IP(I),  I-1,N) 


n o 


1006 


* 

• * 

* 

1007 

* 

* 

1008 

* 

* 

c 

c ***** 


* 


FORMAT  (/'  ROW/COLUMN  ORDERING  FROM  OORV:  '/ 

/'  P IP  ' 

/'  I ROW/COL  ORDERING  INVERSE  ORDERING  ' 
/(I3,  110,  120)) 

1?  (CASE.EQ.2)  PRINT  1007,  (IA(I),  I-1,N),  U{N+1 ) 
FORMAT  (/'  REORDERED  COEFFICIENT  MATRIX:  '/ 

/'  lA  (INDICES  OF  FIRST  ELEMENTS  IN  ROWS)  ' 
/(10I5)) 

IF  (CASE.EQ.2)  PRItll  1008,  (I , JA(I  ) , A(I  ) , I-1,NZA) 
FORMAT  (/'  JA  A ' 

/'  I COLUMN  INDICES  MATRIX  ' 

(13,  110,  F16.5)) 

Call  SDRV  A* 

PATH  - 1 
CALL  SDRV 

(N,  P,IP,  IA,JA,A,  B,  Z,  NSP,ISP,RSP,  PATH,  FLAG) 
IF  (FLAG.NE.O)  GO  TO  102 


AAAAA  Cdlcul<^t6  efTOr  A AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAA 

SUM  - 0 
DO  4 I-1,N 

4 SUM  » SUM  + ((Z(I)-I)/I)**2 
RMS  = SQRT(SUM/N) 

C 

CAAAAAA  Output  sotutton  atid  error  measure  aaaaaaaaaaaaaaaaaaaaaaaaaaaa 
PRINT  1009,  (Z(I),I*1,N) 

1009  FORMAT  (/'  SOLUTION  FROM  SDRV:  ' 

* /(5F10.5)) 

IF  (RMS.LE.EPS)  PRINT  lOlO,  RMS 

1010  FORMAT  (/'  SOLUTION  CORRECT:  RMS  ERROR  = 1PE8.2) 

IF  (RMS.GT.EPS)  PRINT  1011,  RMS 

1011  FORMAT  (/'  SOLUTION  INCORRECT;  RMS  ERROR  = ',  1PE8. 2) 

C 

5 CONTINUE 
STOP 

C 

CAAAAAA  Error  messages  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 


101 

PRINT  1012,  FLAG 

1012 

FORMAT  (/'  ERROR 
STOP 

IN 

OORV: 

FLAG  = ',  15) 

C 

102 

PRINT  1013,  FLAG 

1013 

FORMAT  (/'  ERROR 
STOP 

END 

IN 

SDRV: 

FLAG  - ',  15) 

! 


PT 


Appendix  4 


Sample  Output  From  Test  Driver 


***  FIVE-POINT  OPERATOR  ON  3 BY  3 GRID 
(ALL  ENTRIES  OF  MATRIX  STORED) 

COEFFICIENT  MATRIX: 

lA  (INDICES  OF  FIRST  ELEMENTS  IN  ROWS) 


1 4 8 

11 

15  20 

JA 

A 

I 

COLUMN  INDICES 

MATRIX 

1 

1 

4.00000 

2 

2 

-1 . 00000 

3 

4 

-1 . 00000 

4 

1 

-1.00000 

5 

2 

4.00000 

6 

3 

-1.00000 

7 

5 

-1.00000 

8 

2 

-1.00000 

9 

3 

4.00000 

10 

6 

-1.00000 

11 

1 

-1.00000 

12 

4 

4.00000 

13 

5 

-1 . 00000 

14 

7 

-1.00000 

15 

2 

-1.00000 

16 

4 

-1.00000 

17 

5 

4.00000 

18 

6 

-1.00000 

19 

8 

-1.00000 

20 

3 

-1.00000 

21 

5 

-1 . 00000 

22 

6 

4.00000 

23 

9 

-1.00000 

24 

4 

-1.00000 

25 

7 

4.00000 

?6 

8 

-1.00000 

27 

5 

-1.00000 

28 

7 

-1.00000 

29 

8 

4.00000 

30 

9 

-1.00000 

31 

6 

-1.00000 

32 

8 

-1 . 00000 

33 

9 

4.00000 

RIGHT  HAND  SIDE  B: 

-2.00000  -1. 00000  4.00000  3.00000  0.00000 

7.00000  16.00000  11.00000  22.00000 


ROW/COLUMN  ORDERING  FROM  ODRV 


P IP 

I ROW/COL  ORDERING  INVERSE  ORDERING 


1 

1 

1 

2 

3 

7 

3 

7 

2 

4 

9 

8 

5 

6 

6 

6 

5 

5 

7 

2 

3 

8 

4 

9 

9 

8 

4 

SOLUTION  FROM  SDRV: 

1.00000 

2.00000 

3.00000 

4.00000 

6.00000 

7.00000 

8.00000 

9.00000 

SOLUTION  CORRECT:  RMS  ERROR  = 

1. 39E-08 

***  FIVE-POINT  OPERATOR  ON  3 BY  3 GRID 

(ONLY  ENTRIES  OF 

UPPER  TRUNCLE  STORED) 

COEFFICIENT  MATRIX: 

lA  (INDICES  OF  FIRST  ELEMENTS  IN  ROWS) 

1 4 

7 9 

12  15 

17  19 

JA 

A 

I COLUMN  INDICES 

MATRIX 

1 

1 

4.00000 

2 

2 

-1.00000 

3 

4 

-1.00000 

4 

2 

4.00000 

5 

3 

-1 . 00000 

6 

5 

-1 . 00000 

7 

3 

4.00000 

8 

6 

-1.00000 

9 

4 

4.00000 

10 

5 

-1.00000 

11 

7 

-1.00000 

12 

5 

4.00000 

13 

6 

-1 . 00000 

14 

8 

-1.00000 

15 

6 

4.00000 

16 

9 

-1.00000 

17 

7 

4.00000 

18 

8 

-1.00000 

19 

8 

4.00000 

20 

9 

-1.00000 

21 

9 

4.00000 

RIGHT  HAND 

SIDE  B: 

-2.00000 

-1 . 00000 

4.00000 

3.00000 

7. 00000 

16.00000 

11.00000 

22.00000 

.00000 


1 22 


.00000 


ROW/COLUMN  ORDERING 

FROM  ODRV: 

P 

IP 

I 

ROW/COL  ORDERING 

INVERSE  ORDERING 

1 

1 

1 

2 

3 

7 

3 

7 

2 

4 

9 

8 

5 

6 

6 

6 

5 

5 

7 

2 

3 

8 

4 

9 

9 

8 

4 

REORDERED  COEFFICIENT  MATRU; 

lA 

(INDICES  OF 

FIRST 

ELEMENTS 

IN  ROWS) 

1 4 5 

8 

9 13 

15  18 

JA 

A 

I 

COLUMN  INDICES 

MATRIX 

1 

1 

4.00000 

2 

2 

-1 . 00000 

3 

4 

-1.00000 

4 

2 

4.00000 

5 

2 

-1.00000 

6 

3 

4.00000 

7 

6 

-1.00000 

8 

4 

4.00000 

9 

2 

-1.00000 

10 

4 

-1.00000 

11 

5 

4.00000 

12 

8 

-1.00000 

13 

5 

-1 . 00000 

14 

6 

4.00000 

15 

4 

-1.00000 

16 

7 

4,00000 

17 

8 

-1.00000 

18 

8 

4.00000 

19 

6 

-1.00000 

20 

8 

-1 . 00000 

21 

9 

4.00000 

t 

I SOLUTION  FROM  SDRV: 

1.00000  2.00000  3.00000  4.00000  5.00000 


